home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / SPLINES / EX / PROG7.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  1.8 KB  |  46 lines
  1. // ******************* File: prog7.C *****************
  2. // This program tests the general class TPVolume for
  3. // calculating the spline interpolation of a function
  4. // f(x,y,z,w)=x*y+z*w*w
  5. #include <TPVolume.h>
  6. int main()
  7. {
  8.   // discrete function values on a 4D lattice:
  9.   VecSimple(real) x(7);   // lattice coordinates (in one direction)
  10.   x(1)=-2; x(2)=-1.5; x(3)=-0.5; x(5)=0.5; x(6)=1.5; x(7)=2.0;
  11.   VecSimplest(VecSimple(real)) xs(4); // xs is the 4D lattice coord. info
  12.   for (int i=1; i<=4; i++) { xs(i).redim(7); xs(i)=x; } // same in all dirc.s
  13.  
  14.   Ptv(int) index(4); index.fill(7);   // indexing in 4-dim arrays
  15.   ArrayGenSimple(real) data(index);   // function values
  16.   int j,k,l;
  17.   for (i=1; i<=7; i++)
  18.     for (j=1; j<=7; j++)
  19.       for (k=1; k<=7; k++)
  20.         for (l=1; l<=7; l++)
  21.         { index(1)=i; index(2)=j; index(3)=k; index(4)=l;
  22.           data(index)=x(i)*x(j)+x(k)*sqr(x(l)); }   // f(x,y,z,w)=x*y+z*w*w;
  23.  
  24.   // build knot vectors and spline spaces
  25.   KnotVec knots(5); knots.fill(-2,2); knots.regulate(4);
  26.   SplineSpace space(knots,4);  // cubic spline on knots
  27.   VecSimplest(SplineSpace) spaces(4);  // same space in each direction
  28.   for (i=1; i<=4; i++) spaces(i).redim(space);
  29.  
  30.   TPVolume tpv;
  31.   tpv.redim(spaces);
  32.   tpv.interpolation(xs,data);
  33.   s_o<<"Interpolation in 4D is finished.\n";
  34.  
  35.   // evaluate the spline function at (0.66,-0.9,1.95,0.783):
  36.   Ptv(real) p(4); p(1)=0.66; p(2)=-0.9; p(3)=1.95; p(4)=0.783;
  37.   s_o << "spline at ("; p.print(s_o); s_o << ") = " << tpv(p);
  38.   // ...and the derivative d^2/dzdw at the same point:
  39.   Ptv(int) orders(4); orders(1)=0; orders(2)=0; orders(3)=1; orders(4)=1;
  40.   real result=tpv.derivative(orders,p);
  41.   s_o << "\nand d^2/dzdw =  " << result;
  42.   s_o << "\nExact values are " << 0.66*(-0.9)+1.95*sqr(0.783);
  43.   s_o << " and " << 2*0.783 << '\n';
  44.   return 0;
  45. }
  46.